!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!> Write/read data in the Octave text data format, as used by \c load
!! and \c save in Octave
!!
!! \warning According to the Octave documentation, one should not
!! access data files directly but instead use the Octave library,
!! which has a stable and documented API. Nevertheless, this would
!! require compiling Octave and linking its library, which is
!! sometimes not so obvious. Hence, since the data format is very
!! clean and easy to handle, here we write it ourselves. It could be
!! possible that changes to this module are required when a new Octave
!! release appears, anyway so far this has proven to work very well.
!!
!! For all these subroutines, it is assumed that the file is already
!! open and connected to unit \c fu.
!!
!! Character variables require some care, since one has to deal with:
!! <ul>
!!  <li> single and double quoted strings
!!  <li> strings vs. character arrays
!!  <li> escaping and line termination.
!! </ul>
!! Some documentation can be found <a
!! href="https://www.gnu.org/software/octave/doc/interpreter/Strings.html">here</a>.
!! This module writes double quoted strings and reads both formats.
!<----------------------------------------------------------------------
module mod_octave_io

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_octave_io_constructor, &
   mod_octave_io_destructor,  &
   mod_octave_io_initialized, &
   write_octave,   &
   read_octave,    &
   read_octave_al, &
   real_format,    &
   locate_var

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! private members
 character(len=*), parameter :: &
   real_format = 'e24.15e3' ! format: _+0.***************E+***
!   real_format = 'e42.33e3' ! format: _+0.** 33 **E+*** : real(16)

 character(len=*), parameter :: &
   complex_format = '("(",'//real_format//',",",'//real_format//',")")'

! Module variables

 ! public members
 logical, protected ::               &
   mod_octave_io_initialized = .false.

 interface write_octave
   module procedure write_scal_r, write_vect_r, write_mat_r,  &
                    write_3dim_r, write_4dim_r, write_5dim_r, &
                    write_scal_i, write_vect_i, write_mat_i,  &
                    write_3dim_i, write_4dim_i, write_5dim_i, &
                    write_scal_l, write_vect_l, write_mat_l,  &
                    write_3dim_l, write_4dim_l, write_5dim_l, &
                    write_scal_c, write_mult_c,               &
                    write_scal_z, write_vect_z, write_mat_z,  &
                    write_3dim_z, write_4dim_z, write_5dim_z
 end interface

 interface read_octave
   module procedure read_scal_r, read_vect_r, read_mat_r,  &
                    read_3dim_r, read_4dim_r, read_5dim_r, &
                    read_scal_i, read_vect_i, read_mat_i,  &
                    read_3dim_i, read_4dim_i, read_5dim_i, &
                    read_scal_l, read_vect_l, read_mat_l,  &
                    read_3dim_l, read_4dim_l, read_5dim_l, &
                    read_scal_c
 end interface

 interface read_octave_al
   module procedure read_vect_al_r, read_mat_al_r,  &
                    read_3dim_al_r, & !read_4dim_r, read_5dim_r, &
                    read_vect_al_i, read_mat_al_i,  &
                    read_vect_al_c
!                    read_3dim_i, read_4dim_i, read_5dim_i, &
!                    read_scal_l, read_vect_l, read_mat_l,  &
!                    read_3dim_l, read_4dim_l, read_5dim_l
 end interface

 character(len=*), parameter :: &
   this_mod_name = 'mod_octave_io'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_octave_io_constructor()

  character(len=100) message
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_octave_io_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   write(message,'(A,A)') 'Octave output format: ', real_format
   call info(this_sub_name,this_mod_name,message)

   mod_octave_io_initialized = .true.
 end subroutine mod_octave_io_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_octave_io_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_octave_io_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_octave_io_initialized = .false.
 end subroutine mod_octave_io_destructor

!-----------------------------------------------------------------------
 
 subroutine write_scal_r(x,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(in) :: x
  character(len=*), intent(in) :: var_name
 
  character(len=10) :: out_format
  character(len=*), parameter :: &
    this_sub_name = 'write_octave'

  if (var_name .ne. '') then
    write(fu,'(a,a)')    '# name: ',var_name
  endif
  write(fu,'(a)')      '# type: scalar'

  write(out_format,'(a,a,a)') '(',real_format,')'
  write(fu,out_format) x
  
 end subroutine write_scal_r
 
!-----------------------------------------------------------------------

 subroutine read_scal_r(x,var_name,fu,norewind)
  integer, intent(in) :: fu
  character(len=*), intent(in) :: var_name
  real(wp), intent(out) :: x
  logical, intent(in), optional :: norewind
 
  integer :: ierr
  character(len=1000) :: text
  character(len=*), parameter :: &
    this_sub_name = 'read_octave'

  ! the following subroutine positions the file at the '# name' line
  call locate_var(fu,var_name,ierr,norewind)
  if(ierr.eq.0) then
    read(fu,*) text
    read(fu,*) x
  else
    write(text,'(A,A,A)') 'Problem reading variable "',var_name,'".'
    call warning(this_sub_name,this_mod_name,text)
  endif
  
 end subroutine read_scal_r
 
!-----------------------------------------------------------------------
 
 subroutine write_vect_r(V,rc,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(in) :: V(:)
  character(len=*), intent(in) :: rc, var_name

  integer :: nrow, ncol
  character(len=17) :: out_format
  character(len=100) :: message(2)
  character(len=*), parameter :: &
    this_sub_name = 'write_octave'
 
   if(rc.eq.'c') then
     nrow = product(shape(V))
     ncol = 1
   elseif(rc.eq.'r') then
     nrow = 1
     ncol = product(shape(V))
   else
     write(message(1),'(A)') &
       'Variable rc must be either "r" or "c".'
     write(message(2),'(A,A,A)') &
       'Present value rc="',rc,'".'
     call warning(this_sub_name,this_mod_name,message)
     return
   endif
     
   if (var_name .ne. '') then
     write(fu,'(a,a)')    '# name: ',var_name
   endif
   write(fu,'(a)')      '# type: matrix'
   write(fu,'(a,i7)')   '# rows: ',nrow
   write(fu,'(a,i7)')   '# columns: ',ncol

   ! Note: the repeat count must be strictly positive; on the other
   ! hand it doesn't have to match the data size, so if it's zero we
   ! should set it at least to one
   if(ncol.eq.0) ncol = 1
   write(out_format,'(a,i7,a,a)') '(', ncol, real_format, ')'
   write(fu,out_format) V
  
 end subroutine write_vect_r
 
!-----------------------------------------------------------------------
 
 subroutine read_vect_r(V,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(out) :: V(:)
  character(len=*), intent(in) :: var_name

  integer :: ierr, nrow, ncol
  character(len=100) :: message(3)
  character(len=*), parameter :: &
    this_sub_name = 'read_octave'
 
   call locate_var(fu,var_name,ierr)
   if(ierr.ne.0) then
     write(message(1),'(A,A,A,I3)') &
       'Problems locating "',var_name,'": iostat = ',ierr
     call warning(this_sub_name,this_mod_name,message(1))
   else
     read(fu,'(A)') message(1) ! exclude "# type: matrix"
     read(fu,'(A7,I10)') message(1), nrow ! "# rows:"
     read(fu,'(A10,I10)') message(1), ncol ! "# columns:"
     if((nrow.ne.1).and.(ncol.ne.1)) then
       write(message(1),'(A)') &
         'Variable "',var_name,'" is not a vector.'
       write(message(2),'(A,I10,A,I10)') &
         '  dimension: nrow=',nrow,', ncol=',ncol
       call warning(this_sub_name,this_mod_name,message(1:2))
     elseif(product(shape(V)).ne.(nrow*ncol)) then
       write(message(1),'(A,A,A)') &
         'Dimensions do not agree when reading "',var_name,'":'
       write(message(2),'(A,I10,I10)') &
         'dimension in file: ',ncol,nrow
       write(message(3),'(A,I10)') &
         'internal dimension: ',shape(V)
       call warning(this_sub_name,this_mod_name,message)
     else
       read(fu,*) V
     endif
   endif
  
 end subroutine read_vect_r
 
!-----------------------------------------------------------------------
 
 subroutine read_vect_al_r(V,var_name,fu,norewind)
  integer, intent(in) :: fu
  real(wp), allocatable, intent(out) :: V(:)
  character(len=*), intent(in) :: var_name
  logical, intent(in), optional :: norewind

  integer :: ierr, inrow, incol
  character(len=100) :: message, type
  character(len=*), parameter :: &
    this_sub_name = 'read_octave_al'
 
   call locate_var(fu,var_name,ierr,norewind)
   if(ierr.ne.0) then
     write(message,'(A,A,A,I3)') &
       'Problems locating "',var_name,'": iostat = ',ierr
     call warning(this_sub_name,this_mod_name,message)
   else
     read(fu, '(A7,A)')   message, type  ! "# type:"
     ! A matrix with only one element is saved by octave as a scalar,
     ! which requires a special case.
     if(trim(adjustl(type)).eq."scalar") then
       inrow = 1
       incol = 1
     else
       read(fu, '(A7,I10)') message, inrow ! "# rows:"
       read(fu,'(A10,I10)') message, incol ! "# columns:"
     endif
     allocate( V(inrow*incol) )
     read(fu,*) V
   endif
  
 end subroutine read_vect_al_r
 
!-----------------------------------------------------------------------
 
 subroutine write_mat_r(M,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(in) :: M(:,:)
  character(len=*), intent(in) :: var_name

  integer :: nrow, ncol
  character(len=17) :: out_format
  character(len=*), parameter :: &
    this_sub_name = 'write_octave'
 
   nrow = size(M,1)
   ncol = size(M,2)

   if ( var_name .ne. '') then
     write(fu,'(a,a)')    '# name: ',var_name
   endif
   write(fu,'(a)')      '# type: matrix'
   write(fu,'(a,i7)')   '# rows: ',nrow
   write(fu,'(a,i7)')   '# columns: ',ncol

   ! Note: the repeat count must be strictly positive; on the other
   ! hand it doesn't have to match the data size, so if it's zero we
   ! should set it at least to one
   if(ncol.eq.0) ncol = 1
   write(out_format,'(a,i7,a,a)') '(', ncol, real_format, ')'
   if(ncol.ne.0) & ! otherwise nothing to write
     write(fu,out_format) transpose(M)
  
 end subroutine write_mat_r
 
!-----------------------------------------------------------------------
 
 subroutine read_mat_r(M,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(out) :: M(:,:)
  character(len=*), intent(in) :: var_name

  integer :: ierr, nrow, ncol, i
  character(len=100) :: message(3)
  character(len=*), parameter :: &
    this_sub_name = 'read_octave'
 
   call locate_var(fu,var_name,ierr)
   if(ierr.ne.0) then
     write(message(1),'(A,A,A,I3)') &
       'Problems locating "',var_name,'": iostat = ',ierr
     call warning(this_sub_name,this_mod_name,message(1))
   else
     read(fu,'(A)') message(1) ! exclude "# type: matrix"
     read(fu,'(A7,I10)') message(1), nrow ! "# rows:"
     read(fu,'(A10,I10)') message(1), ncol ! "# columns:"
     if( (size(M,1).ne.nrow).or.(size(M,2).ne.ncol) ) then
       write(message(1),'(A,A,A)') &
         'Dimensions do not agree when reading "',var_name,'":'
       write(message(2),'(A,I10,I10)') &
         'dimension in file: ',ncol,nrow
       write(message(3),'(A,I10,I10)') &
         'internal dimension: ',shape(M)
       call warning(this_sub_name,this_mod_name,message)
     else
       do i=1,nrow
         read(fu,*) M(i,:)
       enddo
     endif
   endif
  
 end subroutine read_mat_r
 
!-----------------------------------------------------------------------
 
 subroutine read_mat_al_r(M,var_name,fu,norewind)
  integer, intent(in) :: fu
  real(wp), allocatable, intent(out) :: M(:,:)
  character(len=*), intent(in) :: var_name
  logical, intent(in), optional :: norewind

  integer :: ierr, inrow, incol, i
  character(len=100) :: message, type
  character(len=*), parameter :: &
    this_sub_name = 'read_octave_al'
 
   call locate_var(fu,var_name,ierr,norewind)
   if(ierr.ne.0) then
     write(message,'(A,A,A,I3)') &
       'Problems locating "',var_name,'": iostat = ',ierr
     call warning(this_sub_name,this_mod_name,message)
   else
     read(fu, '(A7,A)')   message, type  ! "# type:"
     ! A matrix with only one element is saved by octave as a scalar,
     ! which requires a special case.
     if(trim(adjustl(type)).eq."scalar") then
       inrow = 1
       incol = 1
     else
       read(fu, '(A7,I10)') message, inrow ! "# rows:"
       read(fu,'(A10,I10)') message, incol ! "# columns:"
     endif
     allocate( M(inrow,incol) )
     do i=1,inrow
       read(fu,*) M(i,:)
     enddo
   endif
  
 end subroutine read_mat_al_r

!-----------------------------------------------------------------------
 
 subroutine read_3dim_al_r(T,var_name,fu,norewind)
  integer, intent(in) :: fu
  real(wp), allocatable, intent(out) :: T(:,:,:)
  character(len=*), intent(in) :: var_name
  logical, intent(in), optional :: norewind

  integer :: ierr, dims(3)
  character(len=100) :: message, type
  character(len=*), parameter :: &
    this_sub_name = 'read_octave_al'
 
   call locate_var(fu,var_name,ierr,norewind)
   if(ierr.ne.0) then
     write(message,'(A,A,A,I3)') &
       'Problems locating "',var_name,'": iostat = ',ierr
     call warning(this_sub_name,this_mod_name,message)
   else
     read(fu, '(A7,A)')   message, type  ! "# type:"
     ! A matrix with only one element is saved by octave as a scalar,
     ! which requires a special case.
     if(trim(adjustl(type)).eq."scalar") then
       dims(1) = 1
       dims(2) = 1
       dims(3) = 1
     else
       read(fu,*) message ! # ndims: 3
       read(fu,*) dims
     endif
     allocate( T(dims(1),dims(2),dims(3)) )
     read(fu,*) T
   endif
  
 end subroutine read_3dim_al_r
 
!-----------------------------------------------------------------------
 
! Unfortunately, an assumed size array T(*), which would cover all the
! 3D, 4D, xD cases, would not match the generic interface. Thus, it
! would be necessary to call directly the write_multidim subroutine.
! The following is nicer, provided we don't want really large rank
! arrays.
 subroutine write_3dim_r(T,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(in) :: T(:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: dims(3)
  character(len=10) :: out_format
  character(len=*), parameter :: &
    this_sub_name = 'write_octave'

   dims = shape(T)
   write(fu,'(a,a)')     '# name: ',var_name
   write(fu,'(a)')       '# type: matrix'
   write(fu,'(a,i7)')    '# ndims: ',3
   write(fu,'(3i7)')     dims

   write(out_format,'(a,a,a)') '(', real_format, ')'
   write(fu,out_format) T
 
 end subroutine write_3dim_r
 
!-----------------------------------------------------------------------

 subroutine read_3dim_r(T,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(out) :: T(:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: ierr, dims(3)
  character(len=100) :: message(3)
  character(len=*), parameter :: &
    this_sub_name = 'read_octave'

   call locate_var(fu,var_name,ierr)
   if(ierr.ne.0) then
     write(message(1),'(A,A,A,I3)') &
       'Problems locating "',var_name,'": iostat = ',ierr
     call warning(this_sub_name,this_mod_name,message(1))
   else
     read(fu,*) message(1)
     read(fu,*) message(1)
     read(fu,*) dims
     if( (size(T,1).ne.dims(1)).or. &
         (size(T,2).ne.dims(2)).or. &
         (size(T,3).ne.dims(3)) ) then
       write(message(1),'(A,A,A)') &
         'Dimensions do not agree when reading "',var_name,'":'
       write(message(2),'(A,I10,I10,I10)') &
         'dimension in file: ',dims(1),dims(2),dims(3)
       write(message(3),'(A,I10,I10,I10)') &
         'internal dimension: ',shape(T)
       call warning(this_sub_name,this_mod_name,message)
     else
       read(fu,*) T
     endif
   endif
 
 end subroutine read_3dim_r
 
!-----------------------------------------------------------------------

 subroutine write_4dim_r(T,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(in) :: T(:,:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: dims(4)
  character(len=10) :: out_format
  character(len=*), parameter :: &
    this_sub_name = 'write_octave'

   dims = shape(T)
   write(fu,'(a,a)')     '# name: ',var_name
   write(fu,'(a)')       '# type: matrix'
   write(fu,'(a,i7)')    '# ndims: ',4
   write(fu,'(3i7)')     dims

   write(out_format,'(a,a,a)') '(', real_format, ')'
   write(fu,out_format) T
 
 end subroutine write_4dim_r
 
!-----------------------------------------------------------------------

 subroutine read_4dim_r(T,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(out) :: T(:,:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: ierr, dims(4)
  character(len=100) :: message(3)
  character(len=*), parameter :: &
    this_sub_name = 'read_octave'

   call locate_var(fu,var_name,ierr)
   if(ierr.ne.0) then
     write(message(1),'(A,A,A,I3)') &
       'Problems locating "',var_name,'": iostat = ',ierr
     call warning(this_sub_name,this_mod_name,message(1))
   else
     read(fu,*) message(1)
     read(fu,*) message(1)
     read(fu,*) dims
     if( (size(T,1).ne.dims(1)).or. &
         (size(T,2).ne.dims(2)).or. &
         (size(T,3).ne.dims(3)).or. &
         (size(T,4).ne.dims(4)) ) then
       write(message(1),'(A,A,A)') &
         'Dimensions do not agree when reading "',var_name,'":'
       write(message(2),'(A,I10,I10,I10,I10)') &
         'dimension in file: ',dims(1),dims(2),dims(3),dims(4)
       write(message(3),'(A,I10,I10,I10,I10)') &
         'internal dimension: ',shape(T)
       call warning(this_sub_name,this_mod_name,message)
     else
       read(fu,*) T
     endif
   endif
 
 end subroutine read_4dim_r
 
!-----------------------------------------------------------------------

 subroutine write_5dim_r(T,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(in) :: T(:,:,:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: dims(5)
  character(len=10) :: out_format
  character(len=*), parameter :: &
    this_sub_name = 'write_octave'

   dims = shape(T)
   write(fu,'(a,a)')     '# name: ',var_name
   write(fu,'(a)')       '# type: matrix'
   write(fu,'(a,i7)')    '# ndims: ',5
   write(fu,'(3i7)')     dims

   write(out_format,'(a,a,a)') '(', real_format, ')'
   write(fu,out_format) T
 
 end subroutine write_5dim_r
 
!-----------------------------------------------------------------------

 subroutine read_5dim_r(T,var_name,fu)
  integer, intent(in) :: fu
  real(wp), intent(out) :: T(:,:,:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: ierr, dims(5)
  character(len=100) :: message(3)
  character(len=*), parameter :: &
    this_sub_name = 'read_octave'

   call locate_var(fu,var_name,ierr)
   if(ierr.ne.0) then
     write(message(1),'(A,A,A,I3)') &
       'Problems locating "',var_name,'": iostat = ',ierr
     call warning(this_sub_name,this_mod_name,message(1))
   else
     read(fu,*) message(1)
     read(fu,*) message(1)
     read(fu,*) dims
     if( (size(T,1).ne.dims(1)).or. &
         (size(T,2).ne.dims(2)).or. &
         (size(T,3).ne.dims(3)).or. &
         (size(T,4).ne.dims(4)).or. &
         (size(T,5).ne.dims(5)) ) then
       write(message(1),'(A,A,A)') &
         'Dimensions do not agree when reading "',var_name,'":'
       write(message(2),'(A,I10,I10,I10,I10,I10)') &
         'dimension in file: ',dims(1),dims(2),dims(3), &
                               dims(4),dims(5)
       write(message(3),'(A,I10,I10,I10,I10,I10)') &
         'internal dimension: ',shape(T)
       call warning(this_sub_name,this_mod_name,message)
     else
       read(fu,*) T
     endif
   endif
 
 end subroutine read_5dim_r

!-----------------------------------------------------------------------

 subroutine write_scal_i(x,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(in) :: x
  character(len=*), intent(in) :: var_name
 
   call write_scal_r(real(x,wp),var_name,fu)

 end subroutine write_scal_i

!-----------------------------------------------------------------------

 subroutine read_scal_i(x,var_name,fu,norewind)
  integer, intent(in) :: fu
  integer, intent(out) :: x
  character(len=*), intent(in) :: var_name
  logical, intent(in), optional :: norewind
 
  real(wp) :: x_r

   call read_scal_r(x_r,var_name,fu,norewind)
   x = int(x_r)

 end subroutine read_scal_i

!-----------------------------------------------------------------------
 
 subroutine write_vect_i(V,rc,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(in) :: V(:)
  character(len=*), intent(in) :: rc, var_name

   call write_vect_r(real(V,wp),rc,var_name,fu)
  
 end subroutine write_vect_i
 
!-----------------------------------------------------------------------
 
 subroutine read_vect_i(V,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(out) :: V(:)
  character(len=*), intent(in) :: var_name

  real(wp) :: V_r(size(V,1))

   call read_vect_r(V_r,var_name,fu)
   V = int(V_r)
  
 end subroutine read_vect_i
 
!-----------------------------------------------------------------------
 
 subroutine read_vect_al_i(V,var_name,fu,norewind)
  integer, intent(in) :: fu
  integer, allocatable, intent(out) :: V(:)
  character(len=*), intent(in) :: var_name
  logical, intent(in), optional :: norewind

  real(wp), allocatable :: V_r(:)

   call read_vect_al_r(V_r,var_name,fu,norewind)
   allocate( V(size(V_r)) )
   V = int(V_r)
  
 end subroutine read_vect_al_i
 
!-----------------------------------------------------------------------
 
 subroutine write_mat_i(M,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(in) :: M(:,:)
  character(len=*), intent(in) :: var_name

   call write_mat_r(real(M,wp),var_name,fu)

 end subroutine write_mat_i
 
!-----------------------------------------------------------------------
 
 subroutine read_mat_i(M,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(out) :: M(:,:)
  character(len=*), intent(in) :: var_name

  real(wp) :: M_r(size(M,1),size(M,2))

   call read_mat_r(M_r,var_name,fu)
   M = int(M_r)
  
 end subroutine read_mat_i
 
!-----------------------------------------------------------------------

 subroutine read_mat_al_i(M,var_name,fu,norewind)
  integer, intent(in) :: fu
  integer, allocatable, intent(out) :: M(:,:)
  character(len=*), intent(in) :: var_name
  logical, intent(in), optional :: norewind

  real(wp), allocatable :: M_r(:,:)

   call read_mat_al_r(M_r,var_name,fu,norewind)
   allocate( M(size(M_r,1),size(M_r,2)) )
   M = int(M_r)
  
 end subroutine read_mat_al_i
 
!-----------------------------------------------------------------------

 subroutine write_3dim_i(T,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(in) :: T(:,:,:)
  character(len=*), intent(in) :: var_name

   call write_3dim_r(real(T,wp),var_name,fu)

 end subroutine write_3dim_i
 
!-----------------------------------------------------------------------

 subroutine read_3dim_i(T,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(out) :: T(:,:,:)
  character(len=*), intent(in) :: var_name

  real(wp) :: T_r(size(T,1),size(T,2),size(T,3))

   call read_3dim_r(T_r,var_name,fu)
   T = int(T_r)
  
 end subroutine read_3dim_i
 
!-----------------------------------------------------------------------

 subroutine write_4dim_i(T,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(in) :: T(:,:,:,:)
  character(len=*), intent(in) :: var_name

   call write_4dim_r(real(T,wp),var_name,fu)

 end subroutine write_4dim_i
 
!-----------------------------------------------------------------------

 subroutine read_4dim_i(T,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(out) :: T(:,:,:,:)
  character(len=*), intent(in) :: var_name

  real(wp) :: T_r(size(T,1),size(T,2),size(T,3),size(T,4))

   call read_4dim_r(T_r,var_name,fu)
   T = int(T_r)
  
 end subroutine read_4dim_i
 
!-----------------------------------------------------------------------

 subroutine write_5dim_i(T,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(in) :: T(:,:,:,:,:)
  character(len=*), intent(in) :: var_name

   call write_5dim_r(real(T,wp),var_name,fu)

 end subroutine write_5dim_i
 
!-----------------------------------------------------------------------

 subroutine read_5dim_i(T,var_name,fu)
  integer, intent(in) :: fu
  integer, intent(out) :: T(:,:,:,:,:)
  character(len=*), intent(in) :: var_name

  real(wp) :: T_r(size(T,1),size(T,2),size(T,3),size(T,4),size(T,5))

   call read_5dim_r(T_r,var_name,fu)
   T = int(T_r)
  
 end subroutine read_5dim_i
 
!-----------------------------------------------------------------------

 subroutine write_scal_l(x,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(in) :: x
  character(len=*), intent(in) :: var_name
 
  integer :: x_i

   if(x.eqv..true.) then
     x_i = 1
   else
     x_i = 0
   endif
   call write_scal_i(x_i,var_name,fu)

 end subroutine write_scal_l

!-----------------------------------------------------------------------

 subroutine read_scal_l(x,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(out) :: x
  character(len=*), intent(in) :: var_name
 
  integer :: x_i

   call read_scal_i(x_i,var_name,fu)
   if(x_i.eq.1) then
     x = .true.
   else
     x = .false.
   endif

 end subroutine read_scal_l

!-----------------------------------------------------------------------
 
 subroutine write_vect_l(V,rc,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(in) :: V(:)
  character(len=*), intent(in) :: rc,var_name

  integer :: V_i(size(V,1))

   where(V)
     V_i = 1
   elsewhere
     V_i = 0
   endwhere
   call write_vect_i(V_i,rc,var_name,fu)
  
 end subroutine write_vect_l
 
!-----------------------------------------------------------------------
 
 subroutine read_vect_l(V,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(out) :: V(:)
  character(len=*), intent(in) :: var_name

  integer :: V_i(size(V,1))

   call read_vect_i(V_i,var_name,fu)
   where(V_i.eq.1)
     V = .true.
   elsewhere
     V = .false.
   endwhere
  
 end subroutine read_vect_l
 
!-----------------------------------------------------------------------
 
 subroutine write_mat_l(M,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(in) :: M(:,:)
  character(len=*), intent(in) :: var_name

  integer :: M_i(size(M,1),size(M,2))

   where(M)
     M_i = 1
   elsewhere
     M_i = 0
   endwhere
   call write_mat_i(M_i,var_name,fu)
  
 end subroutine write_mat_l
 
!-----------------------------------------------------------------------
 
 subroutine read_mat_l(M,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(out) :: M(:,:)
  character(len=*), intent(in) :: var_name

  integer :: M_i(size(M,1),size(M,2))

   call read_mat_i(M_i,var_name,fu)
   where(M_i.eq.1)
     M = .true.
   elsewhere
     M = .false.
   endwhere
  
 end subroutine read_mat_l
 
!-----------------------------------------------------------------------
 
 subroutine write_3dim_l(T,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(in) :: T(:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: T_i(size(T,1),size(T,2),size(T,3))

   where(T)
     T_i = 1
   elsewhere
     T_i = 0
   endwhere
   call write_3dim_i(T_i,var_name,fu)
  
 end subroutine write_3dim_l
 
!-----------------------------------------------------------------------
 
 subroutine read_3dim_l(T,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(out) :: T(:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: T_i(size(T,1),size(T,2),size(T,3))

   call read_3dim_i(T_i,var_name,fu)
   where(T_i.eq.1)
     T = .true.
   elsewhere
     T = .false.
   endwhere
  
 end subroutine read_3dim_l
 
!-----------------------------------------------------------------------
 
 subroutine write_4dim_l(T,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(in) :: T(:,:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: T_i(size(T,1),size(T,2),size(T,3),size(T,4))

   where(T)
     T_i = 1
   elsewhere
     T_i = 0
   endwhere
   call write_4dim_i(T_i,var_name,fu)
  
 end subroutine write_4dim_l
 
!-----------------------------------------------------------------------
 
 subroutine read_4dim_l(T,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(out) :: T(:,:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: T_i(size(T,1),size(T,2),size(T,3),size(T,4))

   call read_4dim_i(T_i,var_name,fu)
   where(T_i.eq.1)
     T = .true.
   elsewhere
     T = .false.
   endwhere
  
 end subroutine read_4dim_l
 
!-----------------------------------------------------------------------
 
 subroutine write_5dim_l(T,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(in) :: T(:,:,:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: T_i(size(T,1),size(T,2),size(T,3),size(T,4),size(T,5))

   where(T)
     T_i = 1
   elsewhere
     T_i = 0
   endwhere
   call write_5dim_i(T_i,var_name,fu)
  
 end subroutine write_5dim_l
 
!-----------------------------------------------------------------------
 
 subroutine read_5dim_l(T,var_name,fu)
  integer, intent(in) :: fu
  logical, intent(out) :: T(:,:,:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: T_i(size(T,1),size(T,2),size(T,3),size(T,4),size(T,5))

   call read_5dim_i(T_i,var_name,fu)
   where(T_i.eq.1)
     T = .true.
   elsewhere
     T = .false.
   endwhere
  
 end subroutine read_5dim_l
 
!-----------------------------------------------------------------------

 subroutine write_scal_z(x,var_name,fu)
  integer, intent(in) :: fu
  complex(wp), intent(in) :: x
  character(len=*), intent(in) :: var_name
 
  character(len=2+len(complex_format)) :: out_format
  character(len=*), parameter :: &
    this_sub_name = 'write_octave'

   if (var_name .ne. '') then
     write(fu,'(a,a)') '# name: ',var_name
   endif
   write(fu,'(a)')     '# type: complex scalar'

   write(out_format,'(a,a,a)') '(',complex_format,')'
   write(fu,out_format) x
  
 end subroutine write_scal_z

!-----------------------------------------------------------------------
 
 subroutine write_vect_z(V,rc,var_name,fu)
  integer, intent(in) :: fu
  complex(wp), intent(in) :: V(:)
  character(len=*), intent(in) :: rc, var_name

  integer :: nrow, ncol
  character(len=2+7+len(complex_format)) :: out_format
  character(len=100) :: message(2)
  character(len=*), parameter :: &
    this_sub_name = 'write_octave'
 
   if(rc.eq.'c') then
     nrow = product(shape(V))
     ncol = 1
   elseif(rc.eq.'r') then
     nrow = 1
     ncol = product(shape(V))
   else
     write(message(1),'(A)') &
       'Variable rc must be either "r" or "c".'
     write(message(2),'(A,A,A)') &
       'Present value rc="',rc,'".'
     call warning(this_sub_name,this_mod_name,message)
     return
   endif
     
   if (var_name .ne. '') then
     write(fu,'(a,a)')  '# name: ',var_name
   endif
   write(fu,'(a)')      '# type: complex matrix'
   write(fu,'(a,i7)')   '# rows: ',nrow
   write(fu,'(a,i7)')   '# columns: ',ncol

   if(ncol.eq.0) ncol = 1
   write(out_format,'(a,i7,a,a)') '(', ncol, complex_format, ')'
   write(fu,out_format) V
  
 end subroutine write_vect_z
 
!-----------------------------------------------------------------------

 subroutine write_mat_z(M,var_name,fu)
  integer, intent(in) :: fu
  complex(wp), intent(in) :: M(:,:)
  character(len=*), intent(in) :: var_name

  integer :: nrow, ncol
  character(len=2+7+len(complex_format)) :: out_format
  character(len=*), parameter :: &
    this_sub_name = 'write_octave'
 
   nrow = size(M,1)
   ncol = size(M,2)

   if ( var_name .ne. '') then
     write(fu,'(a,a)')  '# name: ',var_name
   endif
   write(fu,'(a)')      '# type: complex matrix'
   write(fu,'(a,i7)')   '# rows: ',nrow
   write(fu,'(a,i7)')   '# columns: ',ncol

   if(ncol.eq.0) ncol = 1
   write(out_format,'(a,i7,a,a)') '(', ncol, complex_format, ')'
   if(ncol.ne.0) & ! otherwise nothing to write
     write(fu,out_format) transpose(M)
  
 end subroutine write_mat_z
 
!-----------------------------------------------------------------------

 subroutine write_3dim_z(T,var_name,fu)
  integer, intent(in) :: fu
  complex(wp), intent(in) :: T(:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: dims(3)
  character(len=2+len(complex_format)) :: out_format
  character(len=*), parameter :: &
    this_sub_name = 'write_octave'

   dims = shape(T)
   write(fu,'(a,a)')     '# name: ',var_name
   write(fu,'(a)')       '# type: complex matrix'
   write(fu,'(a,i7)')    '# ndims: ',3
   write(fu,'(3i7)')     dims

   write(out_format,'(a,a,a)') '(', complex_format, ')'
   write(fu,out_format) T
 
 end subroutine write_3dim_z
 
!-----------------------------------------------------------------------

 subroutine write_4dim_z(T,var_name,fu)
  integer, intent(in) :: fu
  complex(wp), intent(in) :: T(:,:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: dims(4)
  character(len=2+len(complex_format)) :: out_format
  character(len=*), parameter :: &
    this_sub_name = 'write_octave'

   dims = shape(T)
   write(fu,'(a,a)')     '# name: ',var_name
   write(fu,'(a)')       '# type: complex matrix'
   write(fu,'(a,i7)')    '# ndims: ',4
   write(fu,'(4i7)')     dims

   write(out_format,'(a,a,a)') '(', complex_format, ')'
   write(fu,out_format) T
 
 end subroutine write_4dim_z
 
!-----------------------------------------------------------------------

 subroutine write_5dim_z(T,var_name,fu)
  integer, intent(in) :: fu
  complex(wp), intent(in) :: T(:,:,:,:,:)
  character(len=*), intent(in) :: var_name

  integer :: dims(5)
  character(len=2+len(complex_format)) :: out_format
  character(len=*), parameter :: &
    this_sub_name = 'write_octave'

   dims = shape(T)
   write(fu,'(a,a)')     '# name: ',var_name
   write(fu,'(a)')       '# type: complex matrix'
   write(fu,'(a,i7)')    '# ndims: ',5
   write(fu,'(5i7)')     dims

   write(out_format,'(a,a,a)') '(', complex_format, ')'
   write(fu,out_format) T
 
 end subroutine write_5dim_z
 
!-----------------------------------------------------------------------

 subroutine write_scal_c(s,var_name,fu)
  integer, intent(in) :: fu
  character(len=*), intent(in) :: s
  character(len=*), intent(in) :: var_name
 
  write(fu,'(a,a)')    '# name: ',var_name
  write(fu,'(a)')      '# type: string'
  write(fu,'(a)')      '# elements: 1'
  write(fu,'(a,i7)')   '# length: ',len(trim(s))
  write(fu,'(a)')      trim(s)

 end subroutine write_scal_c

!-----------------------------------------------------------------------

 subroutine read_scal_c(s,var_name,fu,norewind)
  integer, intent(in) :: fu
  character(len=*), intent(out) :: s
  character(len=*), intent(in) :: var_name
  logical, intent(in), optional :: norewind

  integer :: i, ierr
  character(len=100) :: message(3)
  character(len=*), parameter :: &
    this_sub_name = 'read_octave'

   call locate_var(fu,var_name,ierr,norewind)
   if(ierr.ne.0) then
     write(message(1),'(A,A,A,I3)') &
       'Problems locating "',var_name,'": iostat = ',ierr
     call warning(this_sub_name,this_mod_name,message(1))
   else
     read(fu,   '(a)'   ) message(1) ! "# type: string"
     read(fu,'(a11,i10)') message(1), i ! "# elements:"
     if(i.ne.1) &
       call error(this_sub_name,this_mod_name, &
            'The input string is not a scalar.')
     read(fu, '(a9,i10)') message(1), i ! "# length:"
     if(i.gt.len(s)) then
       write(message(1),'(a,a,a)') &
         'String "',var_name,'" will be truncated:'
       write(message(2),'(a,i10,a)') &
         '  fortran string: len = ',len(s),','
       write(message(3),'(a,i10,a)') &
         '  octave  string: len = ',i,'.'
       call warning(this_sub_name,this_mod_name,message(1:3))
     endif
     read(fu,'(a)') s
   endif
  
 end subroutine read_scal_c
 
!-----------------------------------------------------------------------

 subroutine write_mult_c(s,var_name,fu)
 ! Write a character array on multiple lines; each line is trimmed.
  integer, intent(in) :: fu
  character(len=*), intent(in) :: s(:)
  character(len=*), intent(in) :: var_name
 
  integer :: lt, i

  write(fu,'(a,a)')    '# name: ',var_name
  write(fu,'(a)')      '# type: string'
  write(fu,'(a)')      '# elements: 1'
  ! determine the total length: first count the newline characters
  lt = max(size(s)-1,0)
  do i=1,size(s)
    lt = lt + len(trim(s(i)))
  enddo
  write(fu,'(a,i7)')   '# length: ',lt
  do i=1,size(s)
    write(fu,'(a)')    trim(s(i))
  enddo

 end subroutine write_mult_c

!-----------------------------------------------------------------------

 !> Read a multiline string, the size of s is the number of lines
 !!
 !! There various formats for a string, here we try to be rather
 !! general and give an error when the format is not recognized.
 subroutine read_vect_al_c(s,var_name,fu,norewind)
  integer, intent(in) :: fu
  character(len=*), allocatable, intent(out) :: s(:)
  character(len=*), intent(in) :: var_name
  logical, intent(in), optional :: norewind

  integer :: ierr, nels, length, ie
  character(len=100) :: message, type
  character(len=len(s)), allocatable :: s_tmp1(:), s_tmp2(:)
  character(len=*), parameter :: &
    this_sub_name = 'read_octave_al'

   call locate_var(fu,var_name,ierr,norewind)
   if(ierr.ne.0) then
     write(message,'(A,A,A,I3)') &
       'Problems locating "',var_name,'": iostat = ',ierr
     call warning(this_sub_name,this_mod_name,message)
   else

     ! Read "# type: <str_type>": there are string, null_string,
     ! sq_string and maybe more; for the known types no special case
     ! should be required.
     read(fu, '(A7,A)')    message, type  ! "# type:"
     if( (trim(type).ne."string"     ) .and. & ! list the known types
         (trim(type).ne."sq_string"  ) .and. &
         (trim(type).ne."null_string") ) &
       call error(this_sub_name,this_mod_name,   &
         'Unknown string type "'//trim(type)//'"')

     ! The next line can be either "# elements:" or "# ndims:",
     ! presently we support only the first format.
     read(fu, '(A11,I10)') message, nels   ! "# elements:"
     if(trim(message).ne."# elements:") &
       call error(this_sub_name,this_mod_name,   &
         'Unknown string format "'//trim(message)//'"')

     ! A string is composed of many elements, each of which can
     ! include multiple lines. Here we separate the various lines,
     ! regardless whether they belong to one or more elements.
     allocate(s(0)) ! initialization
     do ie=1,nels
       nex_el_do: do ! there can be spaces
         read(fu,'(A9,I10)') message, length ! "# length:"
        if(trim(message).eq."# length:") exit nex_el_do
       enddo nex_el_do
       call read_strings(s_tmp1,length,fu)
       ! This would be much simpler using reallocation on assignment
       allocate(s_tmp2(size(s)+size(s_tmp1))); s_tmp2 = (/ s , s_tmp1 /)
       ! compiler bug: see
       ! http://gcc.gnu.org/bugzilla/show_bug.cgi?id=60560
       ierr = size(s_tmp2)
       deallocate(s); allocate(s(ierr)); s = s_tmp2
       deallocate(s_tmp1,s_tmp2)
     enddo

   endif

 contains

  ! The difficulty here is that octave files count the number of
  ! characters in each string, including newline characters and
  ! possible trailing blanks. To make sure that we read exactly the
  ! same number of characters, even when trailing blanks are present,
  ! we use the "size" specifier in the read statement, which requires
  ! advance='no'. Moreover, to make sure that no character are lost in
  ! input, s_tmp is an automatic variable with length equal to the
  ! total number of characters: this is the case when the whole string
  ! is on one line. This requires implementing this part as a separate
  ! subroutine.
  !
  ! Notice that if a string in the input file is too long, it will
  ! still be truncated during the assignment   s = s_tmp   however it
  ! is important not to loose characters during reading, which would
  ! mess up also the following strings.
  subroutine read_strings(s,strlen,fu)
   use iso_fortran_env, only: iostat_eor
   character(len=*), allocatable, intent(out) :: s(:)
   integer, intent(in) :: strlen
   integer, intent(in) :: fu

   integer :: nc, nl, cnt, ios
   character(len=strlen), allocatable :: s_tmp(:)

    ! The strlen characters can be on a single line or they can be
    ! on various lines; the extreme case happens when they are all
    ! newline characters. s_tmp is dimensioned to accomodate all these
    ! cases.
    allocate(s_tmp(strlen)); s_tmp = ''
    nc = 0 ! character counter
    nl = 0 ! line counter
    read_do: do
     if(nc.ge.strlen) exit read_do ! OK also for empty strings
      nl = nl + 1
      read(fu,'(a)',advance='no',size=cnt,iostat=ios) s_tmp(nl)
      ! We should always read until the end of the line, otherwise
      ! there is something wrong either in the file or in this
      ! subroutine.
      if((ios.ne.0).and.(ios.ne.iostat_eor)) &
        call error(this_sub_name,this_mod_name,                       &
          "This should never happen: something wrong either in the "//&
          "input file or in this subroutine.")
      nc = nc + cnt + 1 ! include the final newline
    enddo read_do

    allocate(s(nl)); s = s_tmp(:nl)
    deallocate(s_tmp)

  end subroutine read_strings

 end subroutine read_vect_al_c
 
!-----------------------------------------------------------------------
 
 subroutine locate_var(fu,var_name,ierr,norewind)
 ! Position the file at the first occurrence of string
 ! Use  norewind.eqv..true.  to avoid rewind
  integer, intent(in) :: fu
  character(len=*), intent(in) :: var_name
  integer, intent(out) :: ierr
  logical, intent(in), optional :: norewind

  logical :: found, rew
  character(len=1000) :: record
  character(len=len(var_name)+8) :: record_var_name
  character(len=*), parameter :: &
    this_sub_name = 'locate_var'

   found = .false.
   record_var_name = '# name: ' // var_name! // '\n'

   ! Check whether rewind is required
   rew = .true. ! default is to rewind the file
   if(present(norewind)) rew = .not.norewind
   if(rew) rewind(fu)

   do
    if(found) exit
     read(fu,'(a)',iostat=ierr) record
     if(ierr.ne.0) then ! either error or end of file
       found = .true.
     elseif(record.eq.record_var_name) then
       found = .true.
       ierr = 0
     endif
   enddo

 end subroutine locate_var

!-----------------------------------------------------------------------

end module mod_octave_io

